clear;clc;

load Factors yyyymm;
TBeg = find(yyyymm==201112); TPort = find(yyyymm==201106); 

load CRSPMonthly permno Return Price ShrOut Volume Excd N T;
ME = Price.*ShrOut/1000; Turn = Volume./ShrOut; clear Price ShrOut Volume;
Return11 = NaN(N,T,'single');
for t = 11:T
    Return11(:,t) = prod(1+Return(:,t-10:t),2)-1;
end
Micro = false(N,T);
for t = 1:T
    idx = Excd(:,t)==1;
    x = ME(idx,t); x(~isfinite(x)) = []; x = sort(x);
    Ns = length(x)/10;
    BrkPts = x([1 floor(Ns*(1:10))]);
    idx = ME(:,t)<=BrkPts(3);
    Micro(idx,t) = true;
end
clear t idx x Ns BrkPts;

load BookValue BE_Monthly; BE = BE_Monthly; clear BE_Monthly;
BM = BE./ME;
load CCMAnnualMonthly AssetsGrowth OperatingProfits;
Prof = OperatingProfits./BE; clear OperatingProfits;
Invst = AssetsGrowth; clear AssetsGrowth;

load Data_CRSPMonthly CommonSampleMonthly; % data from 201112 but 201112 has NaN
CommonSample = false(N,T);
CommonSample(:,TBeg:T) = CommonSampleMonthly;
CommonSample(:,TPort:TBeg) = repmat(CommonSampleMonthly(:,2),1,TBeg-TPort+1);
clear CommonSampleMonthly;

NumPort = 10;
temp = struct('AnomalyX',[],'RetPort',[],'WtPort',[]);
VarNames = {'ME', 'BM', 'Prof', 'Invst', 'Return11', 'Return'};
for UseMicro = 0:1
    load Factors yyyymm;
    if UseMicro
        ValidStocks = (Excd==1 | Excd==3) & CommonSample;
    else
        ValidStocks = ~Micro & (Excd==1 | Excd==3) & CommonSample;
    end

    Port = repmat(temp,6,1);

    for vari = 1:4
        Port(vari).AnomalyX = VarNames(vari);
        AnomalyX = eval(VarNames{vari});

        % change signa of some sorting variables so that hi is always expected high return
        if vari==1 || vari==4
            AnomalyX = -AnomalyX;
        end

        RetPort = NaN(T,NumPort,2,'single');
        WtPort = zeros(N,T,NumPort,2,'single');
        for t = TPort:12:T
            X = AnomalyX(:,t);

            idx = isfinite(X) & Excd(:,t)==1 & isfinite(ME(:,t)) & ValidStocks(:,t); % this is the index of stocks to use in calculating breakpoints
            Xsort = unique(sort(X(idx)));  % use unique to do the correct sorts
            NinPort = length(Xsort)/NumPort;
            BrkPts = Xsort([1 floor(NinPort*(1:NumPort))]);
            BrkPts(1) = -Inf; BrkPts(end) = Inf; % set first breakpoint to -Inf; there can be stocks in portfolios that have X1 lower than the lowest X1 from brekpoint stocks

            for n = 1:NumPort
                for s = t+1:min(t+12,T)
                    idx =  X>= BrkPts(n) & X<BrkPts(n+1) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
                    r = Return(idx,s); me = ME(idx,s-1);
                    WtPort(idx,s-1,n,1) = 1/sum(idx); RetPort(s,n,1) = mean(r);
                    WtPort(idx,s-1,n,2) = me/sum(me); RetPort(s,n,2) = sum(r.*me)/sum(me);
                end
            end
        end
        Port(vari).RetPort = RetPort(TBeg:end,:,:);
        Port(vari).WtPort = WtPort(:,TBeg:end,:,:);
    end

    vari = 5;
    Port(vari).AnomalyX = VarNames(vari);
    AnomalyX = eval(VarNames{vari});
    RetPort = NaN(T,NumPort,2,'single');
    WtPort = zeros(N,T,NumPort,2,'single');
    for t = TPort:T-1
        X = AnomalyX(:,t-1); % One extra lag

        idx = isfinite(X) & Excd(:,t)==1 & isfinite(ME(:,t)) & ValidStocks(:,t); % this is the index of stocks to use in calculating breakpoints
        Xsort = unique(sort(X(idx)));  % use unique to do the correct sorts
        NinPort = length(Xsort)/NumPort;
        BrkPts = Xsort([1 floor(NinPort*(1:NumPort))]);
        BrkPts(1) = -Inf; BrkPts(end) = Inf; % set first breakpoint to -Inf; there can be stocks in portfolios that have X1 lower than the lowest X1 from brekpoint stocks

        for n = 1:NumPort
            s = t+1;
            idx =  X>= BrkPts(n) & X<BrkPts(n+1) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
            r = Return(idx,s); me = ME(idx,s-1);
            WtPort(idx,s-1,n,1) = 1/sum(idx); RetPort(s,n,1) = mean(r);
            WtPort(idx,s-1,n,2) = me/sum(me); RetPort(s,n,2) = sum(r.*me)/sum(me);
        end
    end
    Port(vari).RetPort = RetPort(TBeg:end,:,:);
    Port(vari).WtPort = WtPort(:,TBeg:end,:,:);

    vari = 6;
    Port(vari).AnomalyX = VarNames(vari);
    AnomalyX = eval(VarNames{vari});
    AnomalyX = -AnomalyX; % change signs of REVERSAL so that hi is always expected high return
    RetPort = NaN(T,NumPort,2,'single');
    WtPort = zeros(N,T,NumPort,2,'single');
    for t = TPort:T-1
        X = AnomalyX(:,t);

        idx = isfinite(X) & Excd(:,t)==1 & isfinite(ME(:,t)) & ValidStocks(:,t); % this is the index of stocks to use in calculating breakpoints
        Xsort = unique(sort(X(idx)));  % use unique to do the correct sorts
        NinPort = length(Xsort)/NumPort;
        BrkPts = Xsort([1 floor(NinPort*(1:NumPort))]);
        BrkPts(1) = -Inf; BrkPts(end) = Inf; % set first breakpoint to -Inf; there can be stocks in portfolios that have X1 lower than the lowest X1 from brekpoint stocks

        for n = 1:NumPort
            s = t+1;
            idx =  X>= BrkPts(n) & X<BrkPts(n+1) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
            r = Return(idx,s); me = ME(idx,s-1);
            WtPort(idx,s-1,n,1) = 1/sum(idx); RetPort(s,n,1) = mean(r);
            WtPort(idx,s-1,n,2) = me/sum(me); RetPort(s,n,2) = sum(r.*me)/sum(me);
        end
    end
    Port(vari).RetPort = RetPort(TBeg:end,:,:);
    Port(vari).WtPort = WtPort(:,TBeg:end,:,:);

    yyyymm = yyyymm(TBeg:end);

    if UseMicro
        save Portfolios_MicroYes Port permno yyyymm NumPort -v7.3;
    else
        save Portfolios_MicroNo Port permno yyyymm NumPort -v7.3;
    end
end
